#! /usr/bin/python

import sys
import bisect

import pylab
import numpy
import scipy.signal as signal

################### Open the file of interest here #########################
def extract_data(filename):
    infile = open(filename,'r')
    voltage = []
    times = []
    for line in infile:
        words = line.split()
        volt = float(words[0])
        time = float(words[3])
        times.append(time)
        voltage.append(volt)
    infile.close()
    return voltage,times

data_filename = sys.argv[1]
print "reading in %s" % data_filename
sig, times = extract_data(data_filename)
print "done reading in."


###############################################################################

def spectral_inversion(kernel):
    kernel = -kernel
    kernel[len(kernel)/2] += 1.0
    return kernel

def make_fir_filter(sampling_freq, critical_freq, kernel_window, taps, kind,                         **kwargs):
    nyquist_freq = sampling_freq/2
    critical_freq = numpy.array(critical_freq, dtype = numpy.float64)
    normalized_critical_freq = critical_freq/nyquist_freq

    if not taps % 2: #The order must be even for high and bandpass
        taps += 1

    if kind.lower() in ['low','low pass', 'low_pass']:
        kernel = signal.firwin(taps, normalized_critical_freq,
                               window=kernel_window, **kwargs)

    elif kind.lower() in ['high','high pass', 'high_pass']:
        lp_kernel = signal.firwin(taps, normalized_critical_freq,
                                  window = kernel_window, **kwargs)
        kernel = spectral_inversion(lp_kernel)
          
    elif kind.lower() in ['band','band pass', 'band_pass']:
        lp_kernel = signal.firwin(taps, normalized_critical_freq[0],
                                  window = kernel_window, **kwargs)
        hp_kernel = signal.firwin(taps, normalized_critical_freq[1],
                                  window = kernel_window, **kwargs)
        hp_kernel = spectral_inversion(hp_kernel)
        
        bp_kernel = spectral_inversion(lp_kernel + hp_kernel)
        kernel = bp_kernel
    
    return kernel

def fir_filter(sig, sampling_freq, critical_freq, kernel_window = 'hamming',
                taps = 101, kind = 'band', **kwargs):
    """
    Build a filter kernel of type <kind> and apply it to the signal
    Returns the filtered signal.

    Inputs:
        sig          : an n element sequence
        sampling_freq   : rate of data collection (Hz)
        critical_freq   : high and low cutoffs for filtering 
                        -- for bandpass this is a 2 element seq.
        kernel_window   : a string from the list - boxcar, triang, blackman,
                             hamming, bartlett, parzen, bohman, blackmanharris,                              nuttall, barthann
        taps            : the number of taps in the kernel (integer)
        kind            : the kind of filtering to be performed (high,low,band)
        **kwargs        : keywords passed onto scipy.firwin
    Returns:
        filtered signal : an n element seq
    """

    kernel = make_fir_filter(sampling_freq, critical_freq, kernel_window, taps,                                 kind, **kwargs) 


    return numpy.roll(signal.lfilter(kernel, [1], sig), -taps/2+1)

def slice_everything(times, sig, LFP, min, max):
    low_index = 0
    min_index = bisect.bisect_left(times, min)
    max_index = bisect.bisect_left(times, max)
    return times[min_index:max_index], sig[min_index:max_index], LFP[min_index:max_index]

################################## Plot data ###############################
LFP = fir_filter(sig, 10000, [10, 100], kernel_window = 'hamming',
                taps = 2000, kind = 'band')
#downsample the data for exporting to illustrator
ds_factor=20
times=times[::ds_factor]
sig=sig[::ds_factor]
LFP=LFP[::ds_factor]

#pylab.plot(times,sig,'k-',linewidth=1)
stimes, ssig, sLFP = slice_everything(times, sig, LFP, 36000, 42000)

pylab.plot(stimes,sLFP,'k-',linewidth=2)
pylab.ylim(-0.1,0.1)
pylab.show()

